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Abstract 



Smale's a-theory certifies that Newton iterations will converge quadratically to a solution 
of a square system of analytic functions based on the Newton residual and all higher order 
derivatives at the given point. Shub and Smale presented a bound for the higher order 
derivatives of a system of polynomial equations based in part on the degrees of the equations. 
For a given system of polynomial-exponential equations, we consider a related system of 
polynomial-exponential equations and provide a bound on the higher order derivatives of 
this related system. This bound yields a complete algorithm for certifying solutions to 
polynomial-exponential systems, which is implemented in alphaCertif ied. Examples are 
presented to demonstrate this certification algorithm. 
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1 Introduction 

A map / : C" — > C" is called a square system of polynomial-exponential functions if / is 
polynomial in both the variables Xi, . . . ,a;„ and finitely many exponentials of the form e^^' 
where /? € C. That is, there exists a polynomial system P : C""*"™ — )■ C", analytic functions 
gi, . . . ,gm : C — > C, and integers cti, . . . , am € {1, ■ • ■ i such that 



where each Qi satisfies some linear homogeneous partial differential equation (PDE) with complex 
coefficients. In particular, for each i = 1, . . . , m, there exists a positive integer and a linear 
fimction : C^'+i ^ C such that i,{gi, g'^, . . . , g^^) = 0. 
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Consider the square polynomial-exponential system T : C 



^ C"+'" where 



^(xi, . . . , Xfii yi; ■ ■ ■ 7 Um) 



P{X1, ...,Xn,yi,-- ■ ,2/m) 

yi - giixai) 



(1) 



9rn {-^(Tjn ) 



Since the projection map (x, y) ^ x defines a bijection between the solutions of J'{x, y) = and 
f{x) = 0, we will only consider certifying solutions to square systems of polynomial-exponential 
equations of the form F{x,y) — 0. 

For a square system g : C" — >■ C" of analytic functions, a point x € C" is an approximate 
solution of (7 = if Newton iterations applied to x with respect to g quadratically converge 
immediately to a solution of g ~ 0. The certificate returned by our approach that a point is an 
approximation solution of J-" = is an a-theoretic certificate. In short, a-theory, which started 
for systems of analytic equations in [9j, provides a rigorous mathematical foundation for the 
fact that if the Newton residual at the point is small and the higher order derivatives at the 
point are controlled, then the point is an approximate solution. For polynomial systems, by 
exploiting the fact that there are only finitely many nonzero derivatives, Shub and Smale [5] 
provide a bound on all all of the higher order derivatives. For polynomial-exponential systems, 
our approach uses the structure of J- together with the linear functions ii to bound the higher 
order derivatives. 

Systems of polynomial-exponential functions naturally arise in many applications including 
engineering, mathematical physics, and control theory to name a few. On the other hand, such 
functions are typical solutions to systems of linear partial differential equations with constant 
coefficients. Systems, including ubiquitous functions like sin(a;), cos(x), sinh(a;), and cosh(a;), 
can be equivalently reformulated as systems of polynomial-exponential functions, since these 
functions can be expressed as polynomials involving e^^ for suitable /3 e C. Since computing 
all solutions to such systems is often nontrivial, methods for approximating and certifying some 
solutions for general systems is very important, especially in the aforementioned applications. 

In the rest of this section, we introduce the needed concepts from a-theory. Section [2] formu- 
lates the bounds for the higher order derivatives of polynomial-exponential systems and presents 
a certification algorithm for polynomial-exponential systems. In Section |3j we discuss methods 
for generating numerical approximations to solutions of polynomial-exponential systems. Sec- 
tion |4] describes the implementation of the certification algorithm in alphaCertif ied with 
Section [5] demonstrating the algorithms on a collection of examples. 

1.1 Smale's a-theory 

We provide a summary of the elements of a-theory used in the remainder of the article as well as 
in alphaCertif ied. Hence, this section closely follows [H § 1] expect "polynomial" is replaced 
by "analytic." We focus on square systems, which are systems with the same number of variables 
and functions, with more details provided in [5]. 

Let / : C" ^ C" be a system of analytic functions with zeros V(/) = G C" | /(^ = 0} 
and Df{x) be the Jacobian matrix of / at x. For a point x G C", the point Nf{x) is called the 
Newton iteration of f at x where the map Nf : C" — > C" is defined by 




X otherwise. 



X — Df{x) ^f{x) if Df{x) is invertible. 
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For fc G N, let N'^{x) be the fc*'' Newton iteration of / at x, that is, 

N^f{x) ^ NfO---oNf{x). 

k times 

The following defines an approximate solution of / to be a point which converges quadratically 
in the standard Euclidean norm on C" to a point in V(/). 

Definition 1.1 Let / : C" — > C" be an analytic system. A point x G C" is an approximate 
solution of / = with associated solution ^ G V(/) if, for every /c € N, 

||A^)-(^)-eil< 

Clearly, every solution of / = is an approximate solution of / = 0. Additionally, when 
Df(x) is not invertible, then a point x is an approximate solution of / = if and only if 
X e V(/). When Df{x) is invertible, the results of a-theory provide a certificate that x is an 
approximate solution of / = 0. This certificate is based on the constants a{f,x), /3{f,x), and 
7(/, a;) which are defined as 

a{f,x) 
f3{f,x) ^ 

l{f,x) = 

where D^f{x) is the k^^ derivative of / (see Chap. 5]). 

When Df{x) is not invertible, we define I3{f,x) as zero and 7(/, x) as infinity. The constant 
a{f,x) is then the indeterminate form • oo which is defined based on the value of f{x). If 
f{x) = 0, then a{f,x) is defined as zero, otherwise a{f,x) is defined as infinity. 

The following lemma, which is a conclusion of Theorem 2 of [2J Chap. 8], shows that, when 
X is an approximate solution of / = 0, the distance between x and its associated solution can 
be bounded in terms of /?(/, x). Moreover, this bound can be used to produce a certificate that 
two approximate solutions have distinct associated solutions. 

Lemma 1.2 Let f : C" — !■ C" be an analytic system. If x £ C" is an approximate solution of 
f = with associated solution ^, then 

\\x~^\\<2p{f,x). 

Moreover, if Xi,X2 G C" are approximate solutions o/ / = with associated solutions ^i,^2; 
respectively, then ^ ^2 provided that 

\\xi-X2\\>2if3{f,xi) + f3{f,X2)). 
Proof. Both results immediately follow from the triangle inequality. In particular, 

\\x~^\\ < \\x-Nf{x)\\ + \\Nfix)-a</3{f,x)+^-\\x-^\\ 
yields |ja;-^|| < 2/3(/,a;). Additionally, 

\\xi -X2\\< \\xi - ail + iia - 611 + 116 - 2:211 < 2(/3(/, xi) + /?(/, X2)) + Ui- 611 



/3(/,a;).7(/,x), 

\\x-Nfix)\\^\\Dfix)-'fix)\\, and 
Df{x)-^D>^f{x) 



sup 

fe>2 



fc! 



(2) 
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yields that ^ + 6 when \\x^ - > 2(/3(/,a:;i) + /3(/,2;2)). □ 

The following theorem, called an a-theorem, is a version of Theorem 2 of [2, Chap. 8] which 
shows that the value of a{f,x) can be used to produce a certificate that x is an approximate 
solution of / = 0. 

Theorem 1.3 /// : C" — ?> C" is an analytic system and x € C" with 



a{f,x) < ^ — « 0.157671, 

then X is an approximate solution of f — 0. 

The following theorem, called a robust a-theorem and is a version of Theorem 4 and Remark 
6 of [2 Chap. 8], shows that the value of x) and 7(/, x) can be used to produce a certificate 
that X and another point y have the same associated solution. 

Theorem 1.4 Let f : C" — > C" be an analytic system and x G C" with a{f,x) < 0.03. // 
y e C" such that 

then X and y are both approximate solutions of f = with the same associated solution. 

X It 

Let ttk : C" M" be the real projection map defined by 7t-r{x) = — - — where x is the 



conjugate of x. If / is an analytic system such that Nf{x) = Nf{x) for all x such that Df{x) is 
invertible, then Nf defines a real map, i.e., 7Vj(IR") C M". In particular, if x is an approximate 
solution of / = with associated solution 5, then x is also an approximate solution of / = with 
associated solution ^ and /3{f,x) = j3{f,x). The following proposition, which is a summary of 
the approach in [4j § 2.1], can be used to determine if the associated solution of an approximation 
solution is real. 



Proposition 1.5 Let f : <C^ ^ C" be a polynomial system such that Nf{x) = Nf{x) for all 
X € C" such that Df{x) is invertible. Let x g C" is an approximate solution of f = with 
associated solution ^. 

1. //||2:-7rR(a;)|| > 2p{f,x), then ^ ^ M«. 

2. Ifa{f,x) < 0.03 and\\x-TTR(x)\\ < ^^^^^ , t/ien^eM". 

Proof. Since II a; — x| I = 2||a; — 7rR(a;)|| and /?(/, x) = /?(/, x), Item[l]follows by concluding ^ 7^ ^ 
using Lemma 1.2 Item [2] follows from Theorem 1.4 together with the fact that ttm{x) & K" and 
Nf{W') C M". □ 



1.2 Bounding higher order derivatives 

The constant 7(/, x) defined in ^ yields information regarding the higher order derivatives 
of / evaluated at x. Even though, for polynomial systems, ^{f,x) is actually a maximum of 
finitely many values, it is often computationally difficult to compute exactly. However, in the 
polynomial case, it can be bounded above based in part on the degrees of the polynomials [5]. 
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Due to the nature of polynomial-exponential systems, this bound will be used in our algorithm 
presented in Section [2] for certifying solutions to polynomial-exponential systems. 
Let 5 : C" — >■ C be a polynomial of degree d where g{x) = J2\p\<d ^p^'' ^"^^ 



\p\<d 



is the standard unitarily invariant norm on the homogenization of g. For a polynomial system 
/ : C" ^ C" with f{x) = [h{x),..., Ux)V. we have 



ii/iP-Eii/'^ii 



For a point x G C", define \\x\\l = 1 + WxW^ = 1 + YT,=i ■ 

The following is an affinc version of Propositions 1 and 3 from [8j . 

Proposition 1.6 If g : C" — > C is a polynomial of degree d, then, for all x £ C" and k > 1, 

\g{x)\ < \\g\\ ■ \\x\\i and WD^'gix)]] < d ■ {d - 1) ■ ■ ■ {d - k + I) ■ M ■ 

Let k > 2. Lemma 3 of [8 yields 

d-{d~l)---{d-k + l) \ ^ ^ d^^^d-1) ^ d^ 

Additionally, since ||a;||i > 1, we know ||a;||i "'^ > These facts together with Proposi- 

tion |1.6| yield 



d^/^ ■ \\D''gix)\\ 

< ^( d-{d-i)...id-k + i).M.\\x\\t' \^ 



D^g{x) 


k-1 / 

^ ( 


k\ 





fd-{d-l)---{d-k + l) 
di/2 . A;! 



< 



2\\x\ 



which we summarize in the following proposition. 

Proposition 1.7 If g : C" — )• C is a polynomial of degree d, then, for all x G C" and k > 2, 



D^gix) 



kl 



< 



d^ 



2\\x 



Let / 



~^7i \ rnn 



be a polynomial system with degfi = di. Define D — maxd^ and 

^^{f,x) = max{l, ||/|| • \\Df{x)-'A^a)ix)\\} (3) 
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assuming Df{x) is invertible where 



dY'-\\x\\t-' 



(4) 



Since ^{f,x) > 1, < x) for any fc > 2. 

The foUowing version of Proposition 3 of [8, § 1-3] yields an upper bound for 7(/, x). 

Proposition 1.8 Let f : C" — !■ C" be a polynomial system with deg/i = o?,; and D = maxd^ 
For any x G C" such that Df{x) is invertible, 



2-INIIi 



Proof. For A; > 2, we have 
D,f{x)-^D^f{x) 



k\ 



< 



(11/11 •|p/(x)-iA(,)(x)||)-^ 



\f^ 



j3/2 



11/11 -A:! 

1 

2(A;-1)\ 2[fc-i) 



ll/IP l2-||x||i 



< 



2.\\xh ■ 



□ 



2 Certifying solutions 

Since the bound provided in Proposition 1 1 .81 does not apply to a polynomial-exponential system 
we develop a new bound based on the solutions of linear homogeneous partial differential 

equations. With this bound, algorithms for certifying approximate solutions, distinct associated 

solutions, and real associated solutions of 4J apply to 

Consider g{x) — e^^ for some /3 e C. Clearly, for any fc > 0, \g^''Hx)\ — • |5(a;)|. By 

letting B{x) = \g{x)\ and C — max{l, we have 

|5W(x)| <C'=.i?(x). (5) 
The following lemma shows that a similar bound holds in general. 

Lemma 2.1 Let cq, . . . , c^-i G C, £{xo, . . . , Xr) = x^ — X]i=o ^i^iy ^'^d g : C C be an analytic 
function such that ({g,g', ■ ■ ■ ,g^^^) — and r is minimal with such a property. If 

B{x) = max{|5(a;)|, \g' {x)\, . . . , \g'^''~^\x)\} and C = max{l, |co|, . . . , |cr-i|}, 

then, for any a; G C and k > 0, we have 

^ ^' - [ (2 • C)" • r • B{x) ■ C- ifk>r. 
In particular, Ig^'^Hx)] < (2 • C)''-^ ■ r ■ B{x) ■ C = 2^-^ ■ r ■ ■ B{x). 
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Proof. We know g^^^ = J2l^o Cig'^'^\x). For any fc > r, by differentiation, we know 



i=0 

We will now proceed by induction starting at = r. In particular, 

r— 1 r—1 

\9^'\x)\ <Y.\ci\- \g^\x)\ < B{x) • 1 = r • B{x) ■ C. 

i=0 1=0 

For k> r with p = k — r, we have 

1 — 1 / max{r— 1— p,0} r—1 

\g(''){x)\ < Y.\ci\-\g(^+P\x)\<ci ^ \g^'+^Hx)\ + ^ \g(^+P\a 

i=0 y i=0 i=max{0,i — p} 

< C Ir- B{x) + r- B{x)-C ^ 

\ i—r—p 



< r • 



B{x)-c^ i^+c'^Yy^ 



< 2P-r- B{x) ■ CP+^ 
= (2 • C)''-'- ■ r ■ B{x) ■ C. 

The remaining statement follows from the fact that C > 1 and r > 1. □ 

The following lemma will also be used to deduce our bound. 

Lemma 2.2 If 6o > and ai, di,. . . , am, > 1; then 



I m \ 2(fc-i) m 

sup 5o'^'"'^ + 2^*^-'^ E W^O' < ^0 + 2E"?^i 



K>1 



=1 



Proof. Fix fc > 2. Since 2(A; - 1) > 2 and 4(fc - 1) > 2fc, we know af^ > a?'^ and 
^2(fc-i) > £qj. ^ _ rpjjg lemma now follows since 

/ m \ 2(fc-l) / ™ \ 2(fc-l) 

m 

4(fe-l).2(ft-l) 



> ^^('=-^)+22('=-i)^af'=-^)5,' 

m 

> ^^('=-^)+22('=-i)^af^f. 



□ 



Throughout the remainder of this section, we assume that J" : C"+™ — >■ C""'"™ is a polynomial- 
exponential system such that there exists a polynomial system P : C""*"™ — )■ C", analytic func- 
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tions (7i, 



■,971 



I, and integers cri, . . . , (7,-n G {1 

J~ (xi , . . . , Xri ; yi 5 ■ ■ ■ ; Vra) 



. . . , n} such that 

■ • • J Xn I J/l J ■ • • J 2/m) 



Vm 3m (*^cr^ ) 



(6) 



Also, for i = 1, . . . , n, we define di = degP^ and Z? = maxdi. 

We assume that each gi satisfies some nonzero hnear homogeneous PDE with complex co- 
efficients. For each i = 1, . . . ,m, let be the smallest positive integer such that there exists 
a nonzero linear function £i : C"*^^ — > C with li{gi,g[, . . . ,g^^^) = 0. By construction, the 
coefficient of z^. in ^^(zo, zi, . . . , z^i) must be nonzero. Upon rescaling ii, we will assume that 
this coefficient is one, that is, we have 



fi(zo,zi,...,ZrJ 



(7) 



which yields 



'^^j'^Q^ Ci,jgY' ■ We note that the minimal integer with such a property is 



called the order of gi . 

For example, for nonzero A,/i G C, if gi{x) — e^^ , g2{x) — cos(/ix), and 53(2;) — xsin(a;), 
then the order of is 1,2, and 4, respectively. The corresponding differential equations are 



dg 



dx 

with linear functions 



- - A51 = 0, 



/i^(72 — 0, and 



+ ff3 = 



^i(zo,zi) = zi - Azo, £2(zo,2i,-22) = 2:2 +M^2o, and ^3(20, zi, Z2, 23, Z4) = Z4 + 2z2 + zq- 

The bound obtained in Proposition f.8 depends upon fj,{f,x) defined in (|3| for polynomial 
systems. We extend this to polynomial-exponential systems by defining 



{x,y)) = ma.x\ f, 



DF{x,v)- 



A(d)(x,2/)||P|| 



(8) 



assuming that DJ-{x, y) is invertible. The matrix A(^)(a;, y) is the nx n diagonal matrix defined 
in Q and is the m x m identity matrix. We note that ^ reduces to ([3| when m = 0. 
The following theorem yields a bound for j{J-, {x,y)). 

Theorem 2.3 For i = 1, . . . ,m and z £ C, define 

Bi{z) = max{|5i(z)|, . . . , \gl''"~^\z)\} and = max{l, |ci,o|, • ■ • , |ci,n-i|}- 

Then, for any (x, y) G C"^™ such that DJ-{x, y) is invertible, 

( 2^3/2 



Proof. Let M. = 



{x,y)) < K^, {x,y)) 
A(d)(x,2/)||P|| 



DT{x,y)-^D'Tix,y) 



k\ 



mx,v)\\i 

and k > 2. We have 
< \\DT{x,y)-^M\ 



2^Cfmax{f,r, •B,(x^J} 



(9) 



M-^D''T{x,y) 



< ii{F,{x,y)) 



kl 

M-^D^F{x,y) 
fc! 
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By Proposition |1.7| and Lemma [2?T 



M-'^D''F(x,y) 



k\ 



< 



< 



E 

n 

E 



D''P,{x,y) 



d\'^-\\{x,y)rr'-\\P\\-k\ 



\P,\ 



j3/2 



2(fe-l) 



E 



4=1 



IPII 



^3/2 

2||(2:,2;)||i 



m^.v)\\i 

2(fe-l) 



This yields 



sup 

fe>2 



fc! 



< ^(J", (a;,y))sup 

fe>2 



< ^(J", (a:;,y))sup 

fe>2 



^3/2 

2||(a;,y)lli 

^3/2 



2(fe-l) 



2(fe-l) 



22('=-i) ^ (Cf max{l,r, • 

The result now follows from Lemma 12.21 



□ 



Remark 2.4 When m = 0, the bounds provided in Theorem 2.3 and Proposition 1.8 agree. 

The following is an algorithm to certify approximate solutions of = 0. 
Procedure B = CertifySoln(J^, z) 

Input A polynomial-exponential system T : C""*"™ ^n+m ^ point z € C"+™. 

Output A boolean B which is True if z can be certified as an approximate solution of = 0, 
otherwise, False. 



Begin 



1. If J-{z) — 0, return True, otherwise, if DJ-{z) is not invertible, return False. 

2. Set P :— \\DF{z)^^jF{z)\\ and 7 to be the upper bound for "f{F,z) provided in 
Theorem [Ql 



13 — 3\/17 

3. If /3 • 7 < , return Trtte, otherwise return False. 

The algorithms CertifyDistinctSoln and CertifyRealSoln from |3j apply to polynomial- 
exponential systems using the bound provided in Theorem 2.3 The algorithm CertifyDistinct- 
Soln determines if two approximate solutions have distinct associated solutions. The algorithm 
CertifyRealSoln applies to polynomial-exponential systems F such that Njr{W^^™-) c 
and determines if the associated solution to a given approximate solution is real. 

We conclude this section with a refinement of Theorem [23] applied to polynomial-exponential 
systems depending on exp, sin, cos, sinh, and cosh. This refinement uses the following lemma. 
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Lemma 2.5 // Aq, . . . , A„i > and jji, . . . , ii„i > 2, then 

1 

2\ 2(fc-l) 



sup|Af-) + ;^ 



k-1 



k>2 



k\ 



^ m 



i=l 



Proof. Fix k > 2. Since 2{k — 1) > 2 and fii > 2, we know 
follows from 



/^,\2(fc-l) 



> 



. The lemma 



2{k-l) 



> A 







> A, 



2{k-l) 




+ 



^ - 

m 



2(fe-l) 



2(fe-l) 



A 



2{k-l) 



> A 



i=l 

m 2\2(*;-l) 
2(k-l) , ^ Mi A, 



> A 



2(fe-l) 



i=l 



E 



22 



Mi A, 



□ 



Let a, 6, c, e, ft G Z>o, 5i, ejXkiVpy i^q € C, and tTj,T,, ■09 G {1, ... , n}. The following 

considers the following polynomial-exponential system 



g{xi,.. .,Xn,Ui, . . . ,Ua,Vi, . . . ,Vb,Wi, . . . , Wc,2/l, 



,yd,zi, 



,Ze) 



■,Wc, 


yii 


. 1,.. 


,a 


= 1,.. 


.,b 


= 1,.. 


. ,c 


= 1,.. 


. ,e 


= 1,.. 


.,h 


(TN _ 





(10) 



Corollary 2.6 Let Q he defined as in (10) where P : C — !■ C" is a polynomial system with 
N^n + a + b + c + e + h, di = deg Pi and D — max di . For any X, 6 Cz C, define 

A(A, 61) = max{|A|, |A2 exp(A6')/2|}, 

B(A,6I) =max{|A|,|A2sin(A6')/2|,|A2cos(A6l)/2|}, and 

C{\e) = max{|A|, jA^ sinh(A6')/2|, jA^ cosh(A6')/2|}. 

Then, for any X — {x, u, v, w, y, z) € such that DQ{X) is invertible, 

/ 2^3/2 a b c 



^ig,x)<fiig,x) 



i=l 



fc=l 



■^C(?7p,Xp^J + ^C(K5,a:^J I . (11) 

p=l 9=1 / 
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Proof. Let k > 2. The following table lists the bounds on the higher derivatives together with 
associated quantities A and /i used when applying Lemma |2.5[ 



9{x) 


bound for \g^''\x)\ 


A 


A* 




\0''eMOx)\ 


1^1 


max{2, \0cxp{9x)\} 


sm{6x) 
cos{0x) 


|6''=|max{|sin(6'a;)|,|cos(0a;)|} 


1^1 


max{2, |6'sin(fe)|, |6'cos(6'a;)} 


siDh{6x) 
cosh(0a;) 


\e''\ max{| sinh(6'a;)|, | cosh(fe)|} 


1^1 


max{2, |6'sinh(to)|, 16* cosh(6i.T)} 



The result now follows immediately by modifying the proof of Theorem |2.3| incorporating the 



bounds presented in this table together with Lemma 2.5 Based on Lemma 2.5 the functions 
A, B, and C are one-half of the product of the entries in the A and /i columns. □ 



3 Approximating solutions 



In order to certify that a point is an approximate solution of = 0, where is a polynomial- 
exponential system, one needs to first have a candidate point. In some applications, candidate 
points arise naturally from the formulation of the problem. One systematic approach to yield 
candidate points is to replace each analytic function g.^ by a polynomial gf and solve the resulting 
polynomial system, namely 



•^^{^1 ; • ■ ■ 1 ^ni yi ; • ■ ■ 1 ym) 



P{^X\ , . . . , Xn: yi T • ■ • T Un 

yi -9l{Xa, ) 

ym ~ g^^Xom, ) 



(12) 



When the degree of the polynomial approximations are sufhciently large, the numerical solutions 
of = are candidates for being approximate solutions of = 0. In Section [XT] we discuss 
using regeneration [3 to solve T'^ = 0. 

If a numerical solution of = is not an approximate solution of = 0, one can try to 
apply Newton's method for T directly to these points to possibly yield an approximate solution 
of = 0. Another approach is to construct a homotopy between J-^ and J- and numerically 
approximate the endpoint of the path starting with a solution of J-"^ = 0. We note that neither 
method is guaranteed to yield an approximate solution of J-" = 0. 



3.1 Regeneration and polynomial-exponential systems 

Regeneration |3] solves a polynomial system by using solutions to related, but easier to solve, 
polynomial systems. In particular, we will utilize the linear product |13| structure of in (12 1. 



Suppose that g is a univariate polynomial of degree d. The polynomial y — g(x) has a linear 
product structure of 

{x,y,l) X (a;,l) X ••• X (a;, 1) . 



d— 1 times 



This means y—g{x) can be written as a finite sum of polynomials of the form Li {x, y) ■ ■ ■ L^{x, y) 
where 

Li{x,y) — ay + bix + ci and, for z = 2, . . . , d, Li{x,y) = biX + Ci 
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for some a, bi, Ci £ C. 

For i = 1, . . . , m, let = deg and a.i, hi^i, . . . , e C. Similar to the algorithms proposed 
in [3], we note that the following arguments and proposed algorithm depend on the genericity 



of Gi and bij. Define 



Li.i{x,y) ^ aiy + bi,ix + 1 and, for j = 2, . . . , (a;, y) = fe^ + 1. 

Let v — (vi, . . . , Vm) such that I < I'i < ri. Consider the polynomial systems Qi, 



defined by 



: lUli ■ • ■ •] Urn ) 



P{xi, ...,Xn,yi,.. .,ym) 
Ll^vi{Xai,yi) 



(13) 



For 1 = (1,...,1), we first compute the solutions of Qi = 0. We note that in practice, 
Qi is solved by working intrinsically on the linear space defined by Li^^^ (xo-i , yi) = ••• = 
Lm.,u^{xa-^,ym) = 0. Numerical approximations of these solutions can be obtained using stan- 
dard numerical solving methods for square polynomial systems (see [101 I14j ) including, for 
example, polyhedral homotopies |11| or basic regeneration [3]. 

In order to compute the nonsingular isolated solutions of J^^ — 0, we need to compute 
the nonsingular isolated solutions of Qj^ = for all possible ly. By the theory of coefficient- 
parameter homotopies 0, the nonsingular isolated solutions of Qi, = can be obtained by 
using a homotopy from Qi to Qi, starting with the nonsingular isolated solutions of Qi — 0. 
We note that ii i ^ j such that Ui = aj and Vi, Vj > 1, then Qy = has no solutions. 

After solving = for all possible v, we thus have all nonsingular isolated solutions of 



Vixi 



I Xn, yi , 



P{xi,. . .,xn,yi, . • . ,ym) 

I\7=lLl^j{Xa„yi) 



nj" = l Lrn,j{^ 



,2/m) 



= 0. 



(14) 



The final step is to use a homotopy deforming V to starting with the nonsingular isolated 
solutions of T-" = 0. The finite endpoints of this homotopy form a superset of the isolated 
nonsingular solutions of = 0. 



4 Implementation details 

The certification of polynomial-exponential systems is implemented in alphaCertif ied [5^. The 



systems must be of the form in ( 10 ) where the coefficients of P as well as the constant in the 



argument of exp, sin, cos, sinh, and cosh must be rational complex numbers with the bound for 



7 presented in (111. Due to the nature of exponential functions, the computations are performed 
using arbitrary precision floating point arithmetic. Since floating point errors arising from the 
internal computations are not fully controlled, the results of alphaCertif ied for polynomial- 
exponential systems are said to be soft certified. See [H [S] for more details regarding input 
syntax, internal computations, and output. 
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5 Examples 

The following examples used Bertini JV and alphaCertif ied ^ with a 2.4 GHz Opteron 250 
processor running 64-bit Linux with 8 GB of memory. All files for running these examples can 
be found at the website of the first author. 



5.1 A rigid mechanism 

Consider the algebraic kinematics problem fM" of the inverse kinematics of the RR dyad. The 
RR dyad, which is displayed in Figure [l] consists of two legs of fixed length, say ai and 02, 
which are connected by a pin joint. The mechanism is anchored with a pin joint at the point O, 
which we take as the origin. Given a point E = (ei, 62), the problem is compute the angles 9i 
and 62 so that the end of the second leg is at E. That is, we want to solve f {91,62) = where 



f {01,02) = 

The polynomial-exponential system Q 



ai cos(6'i) + 02 cos((?2) — ei 
ai sin(0i) + 02 sin(02) - 62 



of the form ( 10 ) is 



G{0i,e2,yi,y2,y3,y4) 



aiya + a2yA - e.\ 
aiyi + 022/2 - 62 
yi ~ sm{di) 
y2 - sm{d2) 
y^ - cos(6'i) 
yi - cos(6'2) 



Since 9i only appears in / as arguments of the sine and cosine functions, we can compute 
solutions of / = by using the solutions of a related polynomial system. In particular, consider 
the polynomial system g : — !■ obtained by replacing sm{9i) and cos{9i) with and c^, 
respectively, and adding the Pythagorean identities, namely 



5(Sl,S2,Cl,C2) 



aici -I- a2C2 



ei 

a2S2 — 62 

c?-l 



c: 



1 



Given a solution of 5 = 0, solutions of / = are generated using either the arcsin or arccos 
functions. Moreover, it is easy to verify that, for general a^, Ci € C, g = has two solutions and 
thus / = has two 27r-periodic families of solutions. 

Consider the inverse kinematics problem with ai = 3, 02 = 2, and E = (1,3.5). We used 
Bertini to numerically approximate the two solutions of g = 0. For demonstration, consider 
the two digit rational approximations of the solutions 



1 



1 



Xi = —(65,77,76,-64) and X2 — (95, 32, -30, 95) 



100 



100' 



13 



k 









4.94 • 10-^ 


5.26 • 10--^ 


1 


7.46 • 10-" 


6.29 • lO-'* 


2 


1.21 • IQ-i^ 


8.86 • 10-^« 


3 


3.65 • IQ-^^ 


2.01 • 10-^'^ 


4 


3.56 • 10- 


1.10 • 10-^" 


5 


3.56 • 10-"o 


3.41 • 10-"i 


6 


3.50 • 10-'^^^ 


3.21 • 10-^**^ 


7 


3.44 • 10-5S0 


2.90 • 10-^6-* 



Table 1: Newton residuals for Q 



The certified upper bounds for a(g,Xi) computed by alphaCertif ied using exact rational 
arithmetic and rounded to four digits are 0.0736 and 0.0788, respectively. Hence, Xi and X2 
are both approximate solutions of g = 0. Furthermore, alphaCertif ied certified that the 
associated solutions are distinct and real. 

We now consider two corresponding approximations to solutions of ^ = namely 

Zi = (0.711, 2.261, 0.65, 0.77, 0.76, -0.64) and Z2 = (1.874, 0.324, 0.95, 0.32, -0.30, 0.95). 

The upper bounds for a(g,Zi) computed by alphaCertif ied using 96-bit floating point arith- 
metic and rounded to four digits are 0.1265 and 0.1355, respectively. In order to reduce the 
effect of roundoff errors, we also used 1024-bit floating point arithmetic and obtained the same 
four digit value. Hence, alphaCertif ied has soft certified that Yi and I2 are both approximate 
solutions oiQ = 0. Furthermore, alphaCertif ied has soft certified that the associated solutions 
are distinct and real. Table [T] lists the Newton residuals computed by alphaCertif ied using 
4096-bit precision which demonstrates the quadratic convergence of Newton's method. 

By using Euler's formula, we could alternatively use the polynomial-exponential system 
g' -.C^ ^ of the form (fTol) where 



g'{0i,O2,xi,x2,yi,y2) 



aixi + a2X2 - ei + 162 
aivi + a2V2 - ei - ie2 

XiVi - 1 

X2y2 - 1 
j/i - exp(i6'i) 

2/2 - exp(z6'2) 



and i 



-1. Consider the two points 



Wi = (0.711, 2.261, 0.758 - 0.653i, -0.637 - 0.771i, 0.758 + 0.653i, -0.637 + 0.771i) and 
W2 = (1.874, 0.324, -0.299 - 0.954i, 0.948 - 0.318i, -0.299 + 0.954i, 0.948 + 0M8i). 

The upper bounds for a{g' , Wi) computed by alphaCertif ied using both 96-bit and 1024-bit 
floating point arithmetic and rounded to four digits are 0.1492 and 0.1422, respectively. In 
particular, alphaCertif ied soft certified that Wi and W2 are both approximate solutions of 
Q' = with distinct associated solutions. 

Finally, consider the polynomial system obtained by replacing the sine and cosine functions 
in / with a third and second degree truncated Taylor series approximation, respectively, centered 



14 



at the origin, namely 



r 



2) = 



ai(l 



91/2) 
n/6) 



-02(14 
02(^*2 



9i/2)-ei 



01/6) 



62 



The system of equations f^ = has six solutions and yield six solutions of / = upon deforming 
fP to /. These six solutions split into two groups of three based on the values of sm{9i) and 
cos(9i) corresponding to the two families of solutions of / = 0. 



5.2 A compliant mechanism 



In [12] , Su and McCarthy study a polynomial-exponential system modeling a compliant four-bar 
linkage displayed in Figure 4 of [12j . Upon solving a related polynomial system and applying 
Newton's method, they conclude based on the numerical results that a specific compliant four- 
bar linkage has two stable configurations. We will first use alphaCertif ied to certify that their 
numerical approximations of the two stable configurations are indeed approximate solutions. 
Afterwards, we will use the approaches of Section[3]to recompute these two stable configurations. 
The polynomial-exponential system / : — > C''' modeling a compliant four-bar linkage is 



/(a, 6*1, 6*2, 1^1,1^2) = 



R{a){W2 - Wi) + Gi + ncs{9i) - G2 - r2cs{92) 
R{a){W2 - Wi)vi + rics{9i) ~ r2cs{92)v2 
ki{a -a''-9i+ 9\){vi - 1) + /e2(a - 0° - ^2 + Ol){vi 



^2) 



where 



R{a) 



C0s(q;) 
sin(Q:) 



— sin(Q!j 
cos(a) 



and cs{9) 



cos{9) 
siYi{9) 



We note that each of the first two lines in / consists two functions. Additionally, / is not 
algebraic since X, sin(X), and cos{X) all appear in /, where X = a,9i,92- 
The values for the specific linkage under consider are 



Wi = 



-112.632 
-45.053 



VK2 = 



112.632 
-45.053 



, Gi — 



,G2 = 



100 




fci = 29250, k2 = 5824.29, 61? = 1.4486, 9^ = 0.925, and a° 
with numerical approximations for the stable configurations 



Ai = (-0.216933, 1.448567, 0.924966, 0.610174, 1.094669) and 
A2 = (-1.516473, 0.131930, -0.875993, 1.570656, 1.668379). 



,ri =r2 = 250, 



-0.2169. 



The polynomial-exponential system Q 



of the form ( 10 1 is 



Gia,9i,92,i^i,V2,yi, 



i?(yi,y2)(W2 - Wi) -I- Gi -I- rics(i/3, 1/4) - G2 - r2cs{y5,yfi) 

R{yi,y2){W2 - Wi)vi + rics{y3,y4) - r2cs{y5,ye)iy2 
ki{a - a° - 6»i e^){u-, -1) + k2{a ~ a'' - 82 + e^)(i/i - 1^2) 
yi - sin(Q) 
i/2 — cos(q) 
y3 - sm{ei 
t/4 - cos(6'i 
ys - sin(6l2 j 
ye - cos(6'2) 
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F 


bound for 


approximation of 


bound for 


a(F,Bi) 


a{F, B^) 


P{F,B^) 




-f{F,B^), 


liF,B2) 


Q 


1.66 • 10"^ 


4.27 ■ 10-'^ 


8.08 • 10-'' 


1.06 • 10-^ 


2.05 • 10-* 


4.02 • lO'^ 


Q' 


11.9 


42.5 


8.08 • 10-'^ 


1.06- lO"'* 


1.47-10^ 


4.00 • 10^ 



Table 2: Values obtained for Q and G' at Bi and B2 



z 
w 



where 

R{yi,y2) = and cs{w,z) 

\_yi V2 \ 

Let Bi = {Ai^Yi) where 

Yi = (-0.215236,0.976562,0.992539,0.121925,0.798600,0.601862) and 
Y2 = (-0.998525, 0.0542970, 0.131547, 0.991310, -0.768180, 0.640235). 

The upper bounds for a{Q,Bi) computed by alphaCertif ied using both 96-bit and 1024-bit 
floating point arithmetic and rounded to four digits are 0.0166 and 0.0427, respectively. In 
particular, alphaCertif ied has soft certified that Bi and B2 are both approximate solutions 
of C/ = 0. Furthermore, alphaCertif ied has soft certified that the associated solutions are 
distinct and real. 

The formulation of the polynomial-exponential system can have an adverse effect on certifying 
solutions. For example, consider the polynomial-exponential system Q' : C^^ — >■ C"'^^ obtained 
by replacing the 7*^, 9*^, and 11*'' functions of Q with 

yl + yi- 1, yl + yl- 1, and yf + yl- 1. 

Clearly, every solution of ^ = must also be a solution of Q' = 0. Table [2] compares the bounds 
for a and 7 and the value of f3 for Q and Q' at Bi and B2 computed by alphaCertif ied. This 
table shows that the bounds computed for a{G' , Bi) and jiG' , Bi) are three orders of magnitude 
larger than the bounds computed for a{G,Bi) and j{G,Bi). In particular, due to the larger 
bounds, alphaCertif ied is unable to certify that Bi and B2 are approximate solutions of 
G' = 0. If we replace Bi with Ngi{Bi), then alphaCertif ied is able to soft certify that the 
resulting points are approximate solutions oi G' — using both 96-bit and 1024-bit precision. 

We now consider solving a polynomial system obtained by replacing the sine and cosine 
functions with a fifth and fourth degree truncated Taylor series approximation, respectively, 
centered at the origin. Let the polynomial system P : C^^ — > consists of the first five 
functions in G- In particular, P consists of two linear and three quadratic polynomials and thus 



has total degree of the polynomial Qi, defined in ( |13[ ) has total degree 2^ = 8. 

Since we are using fifth and fourth degree polynomial approximations for the sine and cosine 
functions, respectively, we have = 5 if i is odd and r.^ = 4 if « is even. We picked random 
tti, bij G C for J = 1, . . . , 6 and j ~ 1, . . . , and used Bertini to solve each Q^, — 0. In total, 
this produced numerical approximations to 356 nonsingular isolated solutions of T-" = where 



V is defined in ( 14 ) 



The tracking of the 356 paths from V to the polynomial approximation, G"^ , of G produced 
120 points which became the start points for the homotopy deforming \,o Q . This homotopy 
yielded 93 numerical approximations to solutions of 5 = 0. By using both 96-bit and 1024- 
bit floating point arithmetic, alphaCertif ied soft certified that each of these 93 points are 
indeed approximate solutions with distinct associated solutions Moreover, this computation soft 
certified that 65 have real associated solutions, two of which are the two stable configurations 
computed in |12) . 
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